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ABSTRACT 

The purpose of this article is to show that when dynamically cold, dissipationless 
self-gravitating systems collapse, their evolution is a strong function of the symmetry 
in the initial distribution. We explore with a set of pressure-less homogeneous fluids 
the time-evolution of ellipsoidal distributions and map the depth of potential achieved 
during relaxation as function of initial ellipsoid axis ratios. We then perform a series 
of TV-body numerical simulations and contrast their evolution with the fluid solutions. 
We verify an analytic relation between collapse factor C and particle number N in 
spherical symmetry, such that C oc iV 1 / 3 . We sought a similar relation for axisymmetric 
configurations, and found an empirical scaling relation such that C oc N 1 ^ in these 
cases. We then show that when mass distributions do not respect spherical- or axial- 
symmetry, the ensuing gravitational collapse deepens with increasing particle number 
TV but only slowly: 86% of triaxial configurations may collapse by a factor of no more 
than 40 as N — > oo. For N « 10 5 and larger, violent relaxation develops fully under 
the Lin-Mestel-Shu instability such that numerical iV-body solutions now resolve the 
different initial morphologies adequately. 
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1 INTRODUCTION 

Cold, sub-virial distributions of stars undergo a phase of 
gravitational focusing during which binding energy is redis- 
tributed between them (Lynden-Bell's 'violent relaxation' 
process, see Binney & Tremaine 1987 [= BT+87]). The equi- 
libria established through this process show density profiles 
which, when averaged over spherical shells, approach a de 
Vaucouleurs law applicable to elliptical galaxies (van Albada 
1982, McGlynn 1984). This has since motivated studies of 
galaxy and galactic halo formation by some degree of grav- 
itational relaxation (e.g. Henriksen & Widrow 1997, Wein- 
berg 2001). A hands-on approach to this problem, free of 
geometric constraints, consists in integrating the equations 
of motion with iV-body numerical codes. A crucial step when 
applying results from iV-body experiments to actual galax- 
ies and haloes consists in bridging the gap between simu- 
lation particle numbers and the actual number of stars (or 
generally, mass elements) in galaxies, which still differ by 
five orders of magnitude or more in present-day simulations 
of collisionless dynamics (Athanassoula 2000). It is there- 
fore essential to establish the scaling of iV-body results with 
particle number. 

The following example in spherical symmetry brings the 



problem to focus. A star at rest converges to the centre of 
gravity of a free-falling distribution of mass M in an interval 
of time 



iff = 



3tt 



32G <p(a,0)> 



(1) 



where < p(a,0) >= 4-7rM/3a 3 is the mean density inside 
the star's initial radius, a. (The result holds when stars ac- 
crete at the origin, so that shell-crossing is suppressed, see 
Lynden-Bell 1973.) When the density profile is flat initially, 
all stars converge to a point in a free-fall time. The time- 
dependent gravitational potential along a radial orbit is 



$(r,t) = 



GM 
' 2R 



'*!* = 



-™ (2) 



with r < R, the system radius, and double- differencing with 
respect to r at fixed time yields a measure of the tidal field 
at r. We note that the tide is unbound as collapse proceeds 
and R — > 0. In general one would not expect a flat density 
profile on the scales of a galaxy but rather a heterogeneous 
or clumpy matter distribution. Furthermore fragmentation 
modes develop on all scales in homogeneous, cold distribu- 
tions (Aarseth, Lin & Papaloizou 1988, hereafter ALP+88). 
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Bound clumps would survive violent relaxation if their bind- 
ing energy is high (van Albada 1982, Tsuchiya 1998). How- 
ever the above argument suggests that remnant structures 
would be severely affected by the maximal potential depth 
experienced during collapse (here, when R reaches a min- 
imum). What constrains this maximum, and how does it 
scale with particle number? The survival of bound substruc- 
tures, the mixing of orbits and energy exchanges between 
stars all relate to constraints set on (0). 

We explore these questions both analytically, through 
an idealised model of collapsing self-gravitating pressureless 
fluids, and with a set of numerical A-body calculations to 
bring out discreteness effects. We find that while collapse 
simulations in spherical symmetry reproduces theoretical ex- 
pectations of vigorous infall (ALP+88), small departures 
from sphericity lead to much gentler collapse. Specifically, 
we perform A-body simulations of collapsing spheroids cov- 
ering four decades in particle number N (up to 16 millions) 
using an FFT code to integrate the equations of motion. 
We show empirically that during the collapse of spheroidal 
distributions, the maximum gravitational energy achieved 
scales with particle number as oc N 1 ^ 6 and so becomes infi- 
nite as A — > oo. Triaxial initial configurations offer no such 
scaling with particle number. We find from the pressure-free 
fluid calculations that the maximum gravitational energy 
achieved depends sensitively on the initial axis ratios. We 
survey the parameter space of axis ratios on a mesh of 836 
points integrated with high resolution to quantify energy 
maxima. We then use this to constrain the increase in po- 
tential energy of self-gravitating finite- A systems. Thus for 
instance, unless triaxial mass profiles are initially rounder 
than El, their linear size on average will contract by no 
more than a factor O(20). 

We cover the mathematics in Section 2 before giving 
details of our simulations in the section that follows. The 
implications of our study to scenarios of galaxy formation 
are presented in the last section. 



2 ANALYSIS 

Our overall objective is to find out the maximum potential 
depth haloes and galaxies may attain during violent relax- 
ation. For constant initial density, equation (ji|) shows that 
the first caustic occurs at the unique free-fall time tg. This 
would involve the whole of the galaxy, as opposed to a sub- 
set of stars. It is therefore the profile of choice for our study. 
We justify this choice in part below and later with numerical 
modelling (see Section 4). 



2.1 Results in Spherical Symmetry 

We begin by considering in more details the collapse of uni- 
form density spheres. The results of this subsection are those 
of ALP+88. 

The solution for radial infall of uniform spheres takes 
the parametric form 



sin 2rj + 2r\ 



( SGM\ ' 

{7W) x{t ~ to) 



(3) 



with r(t) = r(0) cos 2 rj, rj(t = t ) = and free-fall is com- 
plete when r\ = 7r/2. If we perturb the density profile so 
p : p a + Sp, Sp > (or, = 0) for r < a (or, > a), the free-fall 
time is now a function of position and no singularity forms. 
From ([j]) we write the new collapse time 



t a : t s x (1 - Sp/p + 0[(8p/pf)) = ts- Stff . 



(4) 



Introducing 



»K*ff) = 2 ~ e 



and inserting in (^), we find on truncating the Taylor ex- 
pansion to second order in Sts 



(8GM\ 1/2 A+ 



(5) 



Since the motion is pressure- less, spherical mass shells re- 
expand radially once they have reached the centre. Those 
that originate outside r = a meanwhile continue inwards. 
This spread in arrival time means that the linear dimension 
of the sphere reaches a positive minimum. This minimum is 
known in terms of the original system size and 8tg : 



W(0) 



W(t s ) r(0) 



r {tfi) 2 \ -2 2 /j. n2/3 

= cos r]\tf[) — sin e « e oc (dig) 



(6) 



where C is the collapse factor, and W oc GM 2 jr is the sys- 
tem's gravitational energy. Note that the definition of C ap- 
plies equally to non-spherical systems; the first equality in 
([]) is valid only for spheres. 

So far we have not specified the form of Sp in the region 
< r < a of the initial configuration, only that it be pos- 
itive. If furthermore Sp is a non-monotonic function of po- 
sition, shell-crossing will occur before any shell has reached 
the centre. This would contribute to smooth out fluctuations 
by orbit mixing, but would not affect infall of the system as 
a whole. In the case of a point-mass realisation of a uni- 
form density stellar system, discreteness introduces Poisso- 
nian noise so that Sol p oc 1/yN, with A the total particle 
number. The ratio (0) now scales with the inverse one-third 
power of A, or 



C oc A 1/3 



(7) 



The density perturbations Sp > leading to (^) are Jeans- 
type fragmentation modes of instability: the enhanced grav- 
ity pulls in the matter which condenses faster at the origin, 
as in the classic Jeans condensation of star-formation stud- 
ies. The scaling relation (Q recovers the solution for a cold 
Newtonian fluid in the limit N — ■> 00, for which r(tg) — > 0. 
Systems with small particle numbers experience relatively 
larger density fluctuations, which act as seeds for fragmen- 
tation modes contributing to halt radial infall. Infall stops 
once orbit crossing occurs in the centre, as encapsulated by 
the spread in radial collapse time (^|). 

The above analysis would not apply to initially cuspy 
profiles, since in that case shell-crossing takes place imme- 
diately at the centre. Analysis with shell-crossing is be- 
yond our scope. However if we view a peaked profile as 
a perturbed, uniform-density distribution, where a large- 
amplitude perturbation is necessary to distinguish it from 
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Poissonian noise, the spread in free-fall times, 5tg, would 
then be even larger, so presumably (|^) is minimised for ini- 
tially flat profiles. 

2.2 Non-Spherical Collapse 

2.2.1 Small deviations 

Pressure- less, self-gravitating oblate or prolate structures 
collapse first down the shortest axis (Lynden-Bell 1964, 
Lin, Mestel & Shu 1965). Hence perturbations breaking the 
spherical symmetry while preserving internal density homo- 
geneity grow in time. Writing 

r(t) = r(0)cos 2 77 f + £ (8) 

where r)(t) is known from (^|) and £(r,6,<j>,t) = vector dis- 
placement of the spherical surface with respect to the radial 
direction f. Lagrangian linear analysis shows growth rates 
for these modes (cf. eq. [37] of ALP+88) 

d 2 g = 4tt Gp n - 1 

At 2 ~ 3 (cos 2 ??)" 3 2n + 1* 1 ' 

where n > 1 is the principal number of an angular decom- 
position of the displacement in spherical harmonics (n = 1 
corresponds to a homogeneous radial contraction) . Thus the 
right-hand side in the above equation is positive and |£| be- 
comes larger in time. Since Q is linear in £, angular and 
radial components of each mode (or, value of n) grow in 
time at the same rate. We may solve numerically for |£| as a 
function of time using (^). However an approximate solution 
is found immediately if we note that for collapse in spherical 
symmetry the time-averaged square cosine, 

2 . 1 
<COS V>~ 2 

gives a mean ratio r(£)/r(0) = 1/2 averaged over iff. Substi- 
tuting this in ([]) and writing |£| = £(t) we find on integrat- 
ing 

m = Co exp ( ^8 ^G Po ^± [t - to] J (10) 

Thus high-order n 2> 1 modes grow faster with increasing 
t-t . 

At the time when £(t) ~ r(t), the linear size of the 
system may be compared with (H) in order to determine 
which type of perturbations develop the fastest for a given 
particle number. For this purpose we truncate to third order 
a Taylor expansion with respect to t — t of ( |lo| ) in the 
limit 77 — > 7r/2 (i.e.,t — t a — *• ts)- In the case of Poissonian 
fluctuations, we find from and (Q) that the fragmentation 
modes in spherical symmetry develop faster than surface 
modes ^ when 

N i/s<(« JIEI + ^i^L +1 ) Vi^V 1 (ii) 

~UV2n+l 2 2n+l I \r[0] J 

For a particle realisation of a uniform-density sphere, the 
statistics will be Poissonian. The surface mode should ini- 
tially rise above the noise level to be effective. We therefore 
set 



r(0) s/N ' 
from which we find a critical particle number, 

i/ 6 = ![ fnEL ^n-^_ 

2 V 277 + 1 2 277 + 1 ' v ' 

such that for N > N c fragmentation modes ('dumpiness') 
halt radial infall before V^/V-seed surface modes have de- 
veloped fully and led to pancaking. If we make n > 1 we 
compute the maximum value possible for N c : 

N c w 4.6 6 ~ 9475 . (13) 

Otherwise said, this N c is the largest possible particle num- 
ber for which discreteness effects (noise) may significantly 
distort the flow of a collapsing sphere through a pancaking 
mode. 

The results (|l^) and ( |l3| ) apply to initially small am- 
plitude deviations from spherical symmetry. For sufficiently 
large deviations from sphericity, linear analysis shows that 
small particle number simulations may yet reproduce the 
pancaking collapse of cold fluids (Lynden-Bell 1964, Lin, 
Mestel & Shu 1965). Consider for example an axisymmetric 
displacement of amplitude £ ~ r(0)/4 mapping a sphere to 
a spheroid of aspect ratio 3:4. For this case we compute from 
(0) N ~ 4 3 /(4.6) 3 ~ 1. Thus for axisymmetric cold distri- 
butions, of initial aspect ratio < 3/4, any sensible simulation 
particle number will adequately reproduce the Lin-Mestel- 
Shu flow. The time-evolution of axis-ratio of collapsing tri- 
axial systems with iV = 10 particles performed by Hozumi 
et al. (1996) shows growth of surface modes (pancaking) in 
agreement with ([ll]) . The initial axis ratios of their systems 
were £o/t"(0) ~ 0.01 and 0.005, or three times the Poisson 
noise level for this number of particles. 

We stress that uncorrelated Poisson noise of a uniform 
spherical distribution is not sufficient in itself to lead to ap- 
preciable flattening during infall, for any particle number. 
Thus the scaling (W) is well recovered from simulations with 
as few as N ~ 10 2 particles (see ALP+88, Boily et al. 1999), 
therefore (j^|) does not invalidate the interpretation of previ- 
ous studies of small- iV collapse simulations in terms of one- 
dimensional radial motion (e.g. van Albada 1982, Aguilar & 
Merritt 1990, Cannizzo & Hollister 1992). 

2.2.2 Large deviations: ellipsoidal figures 

An uniform-density ellipsoid collapses down the minor axis 
before major-axis collapse is complete, followed by re- 
expansion when the fluid is also pressure-less. Since all or- 
bits are synchronous in a homogeneous distribution, phase- 
mixing is minimal. Minor-axis cyclic motion continues while 
the ellipsoid collapses down the major axis, until it too re- 
bounds while the minor-axis assumes a finite value. For a 
perfect fluid, such cycles of collapse/expansion may repeat 
themselves without loss of cohesion. For a fluid made up of 
stars however, the stars exchange kinetic energy and phase 
along their orbits, causing damping. In the case where the 
system is initially oblate spheroidal, axes ai = 02 > 03, 
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Table 1. Axes and velocity components for a triaxial pressureless configuration with initial axes (ai, 02,03) = (1, 5, \) for different 
values of the parameter e. The time t max refers to the moment when the configuration is most compact (maximum W). Integration 
ended at i nna i = 4.398 model units. 

£ tmax Wmax Cll dl T T v Cll dl 

H-10 -3 (at t fina i) (at t final ) 



1 x 10" 


-3 


1.934 


6.21 


5 x 10" 


-4 


1.933 


6.31 


1 x 10" 


-4 


1.930 


6.40 


5 x 10" 


-5 


1.931 


6.41 


1 x 10" 


-5 


1.930 


6.42 


5 x 10" 


-6 


1.930 


6.42 


1 x 10" 


-6 


1.930 


6.42 



1.780 


2.291 


0.257 


0.954 


2.291 


0.259 


0.169 


2.291 


0.262 


0.014 


2.291 


0.262 


0.025 


2.291 


0.263 


0.024 


2.291 


0.263 


0.025 


2.291 


0.263 



0.125 0.358 -0.008 

0.169 0.372 -0.044 

0.206 0.377 -0.074 

0.212 0.381 -0.075 

0.213 0.382 -0.076 

0.213 0.383 -0.075 

0.214 0.381 -0.079 



Boily et al. (1999) found in iV-body simulations the time- 
evolution of the aspect ratio after plane crossing to be ap- 
proximately given by 



a 3 (t) _ a 3 (0) ( ai{t) 



ai(t) 01 (0) Voi(0) 



-1/3 



(14) 



i.e. the spheroid becomes rounder as collapse continues and 
the major axis ai(t) — > 0. In practice only a few minor-axis 
oscillations are detected before phase mixing and violent re- 
laxation lead to equilibrium with little or no streaming pat- 
tern. We want to establish for triaxial initial configurations 
a constraint on the collapse factor C in (jfjj) by solving the 
equations of motion for a pressure-less ellipsoidal fluid inte- 
grated over a timescale for complete relaxation suggested by 
JV-body simulations. The motion of a triaxial uniform ellip- 
soid, of axes ai > a,2 > 03 , is governed by a set of harmonic 
equations (BT+87, Table 2.1) 



At x{ 



(15) 



(16) 



where Xi(t) = a,i(t)/ai(6), v(t) = dxi/ At are dimensionless 
functions of the dimensionless time t = t/tg. The coefficients 
Ai(a) are known from potential theory. For instance we have 

, . 02(0)03(0) F(6,k)-E(9,k) 
U ' af(0) fe 2 sin 3 6» 

with similar definitions for A2,A 3 , and 



k(t) 



ai 



1/2 



lit) 



Ol 



(17) 



Here F, E are incomplete elliptical integrals (see BT+87 for 
details). We may identify the most relevant configurations 
by inspecting the gravitational energy. The self-gravitating 
potential energy W is known for uniform ellipsoids from 



W(a,t) 



3 M 2 F(6,k) 



(18) 



5 ai(i) sin 9 

with the definitions ([n]). The energy W in (^) diverges 
when the ellipsoid collapses to a rod (or spindle) which is 
the case when a — > (ai, 0, 0). It remains finite when two axes 



are non-zero, which includes collapse to a disc. Presumably 
the force and tidal fields are maximised when ellipsoids de- 
velop spindles, together with W. To constrain the tidal field, 
it is therefore sufficient to determine when an ellipsoid forms 
a spindle, or, generally, what maximum value W may reach 
during evolution. We were not able to determine these an- 
alytically and have resorted to a numerical integration of 
the equations of motion. We found it useful to introduce the 
parameter r defined as 



a-z — as 
en — C13 



> 



(19) 



Thus axisymmetric prolate spheroids (02 =£13) all have r = 
0, whereas axisymmetric oblate spheroids (ai = 02) have 
t = 1. Triaxial structures assume intermediate values. 



3 SOLUTIONS FOR PRESSURE-LESS 
ELLIPSOIDS. 

3.1 Method & Tests 

For given axes a = (0,1,0,2,0,3) we may integrate (^) and 
( p^| ) subject to the initial conditions Xi(0) = l,Vi = 0. 
Note that the equations are singular when any of the axes 
vanishes. To integrate through such singularities we enforce 
time-symmetry by reversing the flow: Vi — ► — Vi whenever 
Xi < e, with e a free parameter. We set up a fourth-order 
Runge-Kutta integrator (Press et al. 1992) and varied e from 
lCF 3 down to 1 x 1CP 6 without appreciable differences in the 
integrated global quantities such as maximum W and time 
(see Table hi). However details of the fluid configurations 
converged to good accuracy only when e~5x 10~ 5 or less. 
To quantify the quality of orbit-integration, we computed 
both axial lengths and velocity components for the repre- 
sentative case where (ai,a2,a3) = (1,5,5); then initially 
t = j from (|l^). We computed a similar quantity r v from 
the velocity components (di, 02,03) which we evaluated at 
t = tmax, when W/W(0) is maximum. 

The results are listed in Table |l| for various values of e. 
By comparing the runs of individual components (01, di) 
and those of (j,t v ) with decreasing e, we may conclude 
that both geometric and velocity ellipsoids vary little with 
e. This does not hold for individual components, such as 
major axis, 01. This last quantity must reach a minimum a 
few times larger than e to ensure that the dynamics is re- 
solved properly. Drawing from the results in the Table, we 
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Figure 1. Maximum potential energy Wmax achieved by pressure- less homogeneous ellipsoids of axes a\ = 1 > a 2 > 113 collapsing from 
rest. W has been normalised to its initial value. The Wmax surface is projected in the 02 — (13 plane and three contours are shown, 
W/W = 12, 18 and 24. The maximum values displayed have been truncated to 25. 



set e = 2 x 10 5 or 1.5 x 10 5 in all our integrations as a 
minimal condition to accurate integration. 

In conjunction with the value of e, the choice of timestep 
is crucial: we adopted a time-adaptative scheme tailored to 
the instantaneous free-fall time (|l|) at each step. This allowed 
us to resolve in time increases in potential energy by factors 
up to « 340, while keeping errors below the 1% level (though 
not for axisymmetric systems, see below). 

When integrating equations ( EEl ) and (|l6|), care must be 
taken that the indices (1,2,3) are circulated between each 



axis to identify the current major and minor axes properly, 
and allow the correct evaluation of the force field. All in- 
tegrations were done in three dimensions but we found it 
necessary to enforce symmetry in the potential when treat- 
ing initially axisymmetric configurations in order to prevent 
large numerical errors. With enforced symmetry, we com- 
puted correctly the growth of singular axisymmetric poten- 
tials for a collapse factor reaching w 40. As a test, we inte- 
grated through the singularities formed through major-axis 
collapse of axisymmetric spheroids with initial axes (1, |, |) 
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and (1,1, §)• The error in binding energy at the end of inte- 
gration reached 5% and 6% in each case respectively, which 
we take as reference later when sampling the space of triax- 
ial initial configurations. Note that these errors were accrued 
during a single integration through singularity: the slow di- 
vergence of W — » oo as axes vanish makes agreement with 
theory intractable, however the important quantity for these 
cases must be the time t when W diverges (since the spin- 
dle morphology is known) , which could be evaluated to high 
accuracy. 

In order to determine the sensitivity of our integration 
scheme to the symmetry of the system, we integrated a tri- 
axial configuration with axes = (1, 0.99, |), i.e. nearly oblate 
axisymmetric. In this case integration yielded a maximum 
collapse factor of m 26 and total error accrued for the sys- 
tem binding energy of w 0.5% with the same set-up used 
for the strictly axisymmetric case discussed before. We con- 
cluded from this that the integrator resolves relatively small 
departures from axisymmetry, of the order of one percent in 
axial ratios; and that these small departures are sufficient to 
avoid spindles and large energy errors, while still collapsing 
by appreciable factors. 

We have monitored the three function of 

time to determine whether a spindle forms, which, in view 
of our approximations, is the case when 

a — > (oi, e, e) ; ai > £ : definition of a spindle. (20) 

If this occurs, then a self-gravitating body may assume an 
even larger collapse factor C defined in (^) since our scheme 
does not resolve the dynamics below that scale. The iden- 
tification of configurations of high potential energy is made 
difficult due to the slow divergence of W with vanishing mi- 
nor axes 02, 03 (spindle), or all three axes simultaneously 
(spherical case). By contrast, the case when only one axis 
vanishes (by definition not the major axis) is integrable to 
high accuracy. We therefore computed the logarithmic aver- 
aged length I = (aifl2a3) 1 ' 3 , which we used together with 
( P0|) to ensure that the maximal potential energy computed 
occurred at the same time as the minimum of I. As a final 
precaution, we rejected any integrated solution that accu- 
mulated errors in binding energy exceeding 10%, in view of 
our tests with axisymmetric configurations. When ( |2o| ) does 
not occur, the collapse factor reaches a finite maximum: this 
maximum is then the limit any TV-body collapse calculation 
for this initial geometry may reach, since discreteness effects 
will only increase the growth of kinetic energy and slow down 
collapse. 

3.2 Results 

Anticipating the results from TV-body calculations of Sec- 
tion |E| we integrated (|l^) and (^6|) up to t — tg (t = 1), 
corresponding to one free-fall time (jjj) in spherical symme- 
try. We then explored the parameter space of (02(0), 03(0)) 
by sampling the parameter r uniformly in the interval [0,1], 
for a total of 900 pairs (02,03), while fixing the major axis 
to an initial value ai(0) = 1. We evaluated ( plj| ) and kept the 
largest value, W m ax, found during integration. The results 
are displayed on Fig. |l[ The bottom panel graphs the surface 
of maximum energy for all pairs (02(0), 03(0)). It is striking 
that large islands exist where Ts/Lscx(W) ~ O(10), whereas 



Figure 2. Integrated distribution (in %) for ellipsoidal pressure- 
less uniform fluids as function of the maximum gravitational en- 
ergy achieved during infall. For each 10% interval, we give the 
mean shape parameter <t> (eq. Jl9| ) and its standard deviation 
evaluated from the initial conditions. 

all axisymmetric configurations with r = 1 or (i.e., 02 = 1 
and 02 = 03, respectively) must develop spindles and infinite 
W. The vertical axis has been capped to W m ax /W (0) = 25 
for clarity. Larger values were not prohibited in the course 
of integrating numerically. The presence of fragmented re- 
gions with large W in the plane (02, 03) are indicative of the 
formation of spindles or very compact configurations, in the 
sense of our equation (^) . We note the presence of islands of 
as few as a single point were W m ax/W(Q) > 25, suggesting 
a complex topology. Details of the topology of the energy 
surface W are not important to the main argument and will 
not be pursued further. 

Cases where spindles formed from initially triaxial con- 
figurations turn out to be exceptional. For these cases, a 
repeat of the integration with a smaller e lead to larger 
Max(W) at the time the spindle formed, requiring care- 
ful step-wise integration over a small time interval. The re- 
mainder of the integrated solution, however, was left largely 
unchanged, giving confidence that singularities were cor- 
rectly identified and cured. As stated earlier, we rejected 
all runs that accumulated errors in binding energy through 
full time-integration exceeding 10%: in total some 62 (i.e. 
7.1%) of all initial configurations were rejected on this ba- 
sis. A graph of 5E/E vs Max(W / /M / (0)) showed no clear 
trend with increasing Max(W/W / (0)) for these 62 cases, as 
would have been the case if a single integration through sin- 
gularity at high density was accountable for the bulk of the 
error budget: instead, large numerical errors develop owing 
to repeated integration through singularities, which may oc- 
cur in rapid succession if the initial configuration is signifi- 
cantly non-axisymmetric. Indeed of the 62 cases with signif- 
icant energy errors, 37 initially had axis ratios 03/01 = 1/10 
or lower; two more showed near sphericity, with 02/01 and 
03/ai > 0.999. The remaining 23, however, showed no pecu- 
liarities in their initial values, or in the maximum W com- 
puted. They were, nevertheless, left out of the analysis. 

From Fig. El we may quantify the fractional area in 
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Figure 3. Top panel: scatter plot showing the maximum W 
achieved as function of the initial shape parameter to. Lower 
panel: comparison between the initial morphology (measured by 
t, cf. eq. h9() and the morphology when W is maximal. The inte- 
grated distribution of r is shown in each case. This demonstrates 
a drift towards lower r during infall. 



(02,03) leading to large potential energy and collapse fac- 
tor. The inset gives the projected isocontours in the param- 
eter space. If we assume fair sampling, then the total area 
covered by the highest-level contours is estimated easily us- 
ing a rectangular grid to cover the contoured area. In this 
way we compute a net fraction of w 30% of the total area 
exceeding an increase in W/W(0) of 25. Note that practi- 
cally all triaxial configurations with 02 ^ 0.4 and 03 < 0.2 
likely develop large W . These considerations are quantified 
more accurately by sorting Wmax/W(0) in increasing or- 
der for all pairs (02,03). This gives an integrated distri- 
bution of collapse factor C. Figure ^ plots the integrated 
fraction of 836 solutions as function of Wmax/W(0). Two- 
thirds of these solutions reached a collapse factor C < 22.4, 
and 86% have C < 40. We sought a correlation between 



the morphological parameter r and the maximum potential 
reached. On Fig. |^ we also plot the mean r computed at 
t = for each ten-percentile interval, in increasing order 
of Wmax/W(0). We find a broad trend such that ellipsoidal 
initial conditions with larger r tend to collapse to deeper 
potentials. We may identify for the first of these bins, which 
has < t >= 0.24, the broad \ow-W m ax valley seen on Fig. |l| 
The non-monotonic relation of <r> with W m ax/W(0) may 
be guessed at if we look at a scatter plot of this quantity ver- 
sus r directly, as shown on the top panel of Fig. |3[ The points 
clustered to the bottom left corner of the graph clearly re- 
flect the trend of small r to yield small W m ax/W(0). The 
broad trend we compute for <r> can be guessed from shift- 
ing a horizontal ruler vertically up the W axis: there is a 
suggestion of a gap in the data which account for the dip 
in < r > in the 50-60% interval (when W max /W(0) ^ 20). 
Our conclusions concerning the significance of this gap must 
be moderated by the large deviations about the mean val- 
ues. A more robust signature of dynamical evolution is a 
shift of the distribution of r towards lower values during 
infall. This may be measured by computing the shape pa- 
rameter t at the time when W is maximum, and compared 
to the initial value, To, for a given configuration. The re- 
sult is shown on the bottom panel of Fig. ^. At maximum 
W/W(0), the mean r = 0.29 approximately, down from 0.50 
for the initially uniform distribution. By the definition (|l9|), 
this implies a roundening of the two minor axes, in the same 
fashion as occurs for axisymmetric spheroids, that is 03/02 
increases during infall (see Boily et al. 1999 for a discussion 
of this issue) . If the systems were allowed to virialise, as will 
happen in an iV-body calculation allowing orbit mixing and 
violent relaxation, in contrast to the fluid model presented 
here, the effect measured during infall of the increasing ra- 
tio 03/02, so enhancing axisymmetry, would be offset by the 
onset of radial orbit instability, which is known to enhance 
triaxiality (see e.g. Aguilar & Merritt 1990). We have not 
addressed here the question of which effect would prevail in 
determining the equilibrium of haloes or galaxies. This issue 
would require more realistic density profiles than used here 
and is beyond the scope of the present paper. 

3.3 Summary 

We sum up the results for pressure-less self-gravitating el- 
lipsoidal collapse as follows: 1) The bulk of initially triaxial 
figures does not reach a collapse factor C = W/W(0) ex- 
ceeding 40. We find that 86% of the 836 configurations ex- 
amined collapsed to smaller values; 2) Triaxial figures that 
collapse by small factors (small increase of W/W(0)) tend to 
have preferentially small values of the parameter r (i.e., are 
more axisymmetric, see Fig. ^); 3) the axis ratio 03 to 02 
increases during infall, leading to more axisymmetric con- 
figuration at maximum collapse (Fig. ^); 4) since no orbit 
mixing or fragmentation mode develop in the fluid solutions, 
which would contribute to boost velocity dispersion and 
stop infall , we deduce that the collapse factors C obtained 
are absolute maxima. 

In order to apply correctly these results to galaxies and 
haloes, we must first quantify the impact of discreteness ef- 
fects of finite- N systems as discussed in Section 2. We pro- 
ceed empirically with iV-body numerical calculations. We 
turn first to the task of reproducing the theoretical expecta- 
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tion (pi) for spherically symmetric systems. This is followed 
by a set of calculations of non-spherical initial configura- 
tions, from which we seek trends with particle number to 
compare with the results for perfect fluids of Fig. [U 



4 TV-BODY CALCULATIONS: SET-UP AND 
TESTS 

4.1 Numerical code and units 

Our intention is to cover as wide a range in particle number 
as possible to seek out correlations applicable to galaxies and 
haloes. The nested-grid code SUPERBOX was used (Fellhauer 
et al. 2000). This is an FFT Poisson solver with Cartesian 
grids, of uniform resolution for each cell of a given grid. 
There are three grid levels, each within one another, thus 
enabling enhanced resolution where it is required. This is 
a crucial feature for the problem at hand, since the struc- 
tures collapse by large factors and the density increases ac- 
cordingly. For computational purposes the units were chosen 
such that the total mass and initial radius of the system 
I . however the gravitational constant G — 2. The free-fall 
time (^) for uniform spheres is therefore 

■ = - ~ 0.785398... . (21) 
4 

A version of the code was set up where the integration 
timestep of a leap-frog scheme scales down with the inverse 
square root of the local density maximum during the sim- 
ulation, in order that the timestep St remains in the same 
ratio to the instantaneous dynamical time (or, ts evaluated 
from |l] but with <p> now the density at time t) . The three 
levels of resolution allowed by the code were set such that 
when the system reaches a minimum size, the number of 
particles per cell is of the order of a few, hence interactions 
between particles are resolved approximately on the smallest 
scales. We found in practice that grid resolution smoothes 
out forces between particles and hence leads to a less deep 
potential minimum at the centre of gravity. As we show be- 
low, the net effect of gridding can be easily brought to low 
error levels. 

4.2 Tests in spherical symmetry 

We carried out checks of our set-up in spherical symmetry 
before conducting a survey of particle number and geometry. 
All simulations of violent relaxation proceed from zero ve- 
locity initially. Our results are summarised in Table ^| The 
analysis of Section 2.1 suggests a clear relation (Kepler's) 
between system radius and the time interval Stg before col- 
lapse. From (^) and (|^), 

r(tg - i) « r(0) sin 2 e oc (t B - t) 2/3 . 

We recover this relation for large particle number simula- 
tions shown on Fig. ^. This graphs the evolution of four con- 
stant mass (Lagrangian) shells for a uniform density sphere 
of one million particles. The initial radius r(0) = 1. Be- 
cause the free-fall time is independent of position, the La- 
grangian radii must remain in the same relative ratio to 
one-another; each must converge to the Keplerian regime 



near full collapse. The time of collapse (t — tg = 0) is 
off-scale on the right-hand side of the logarithmic abscis- 
sae on the figure. Two set-ups are illustrated, of low- (left- 
hand) and high-resolution (right-hand) grids. The linear 
high-resolution achieved, I = grid size / number of cells 
= 0.05/64 = 0.0008, is still large when compared with the 
mean inter-particle distance U n t expected from (g), 

Volume at bounce = TV (li„t/2) 3 — 

Initial volume -4- N = r(0) 3 /TV, 

or Unt r; 2 x TV -2 / 3 = 2 x 10~ 4 . If we count the average 
number of particles in a (cubic) cell at the bounce, we find 
~ (l/'ihnt) 3 ~ O(l) particles. The mass distribution is there- 
fore well sampled, and as a result both the constant ratios 
of Lagrangian radii and their match of the Keplerian rela- 
tion are reasonably well recovered. By comparison, for runs 
with reduced number of mesh points, from 128 to 32, we find 
« 100 particles in each cell at the bounce. Poor resolution 
of the mass distribution leads to artificial deviations from 
the Keplerian tracks (see left-hand panels on Fig. ^). The 
Lagrangian radii spread out, which results in shell crossing 
at the centre while the outer shells continue to fall in: this 
causes artificial orbit mixing and a gentler collapse. 

Our computational strategy must therefore ensure that 
the mass profile is well resolved at all times. Because the 
one-dimensional spherical collapse provides the strictest nu- 
merical test of our numerical set-up, we first recover the 
scaling (|?j) for large particle number to refine the code's grid 
resolution. We do this for values of TV ranging from 10 4 to 
16 x 10 6 . The mass distribution is mapped accurately by La- 
grangian radii sorted on concentric shells. We have measured 
the collapse factor (^|) using both the ratio of gravitational 
radius, r g , and two shells enclosing 30% and 60% of the to- 
tal mass. As can be seen from Table [| the collapse factor 
found for given grid resolution and particle number varies 
considerably depending on the choice of Lagrangian radii, as 
well as with the ratio of gravitational energies. However the 
free-fall time is recovered to 1.6% or better. These different 
values obtained for the collapse factor were used to define 
error bars on the averaged quantities. The trend with par- 
ticle number is displayed on Fig. |E|(a). Results in spherical 
symmetry are plotted as circles on the figure: the filled cir- 
cles represent results from this paper, while open circles are 
taken from ALP+88. Each point is the average of data listed 
in Table ^ for given TV, but excluding the low-resolution runs. 
The scaling relation ([?]) shown as the solid line on the figure 
provides a good fit to all data points. 



5 RESULTS FOR TV-BODY SPHEROIDAL 
COLLAPSE (r = OR 1) 

Our results for pressure-less fluids show that initially ax- 
isymmetric distributions (with r = or 1) develop spindles 
and divergent gravitational energy as they collapse. These 
configurations may as a result show significant dependencies 
on particle number in an TV-body realisation of the solution. 
We therefore explore the case of the collapse of spheroidal 
distributions first. 
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Figure 4. Orbits of four Lagrangian mass shells during the infall of a uniform-density sphere. The shells enclose 10, 30, 50 and 80 % of 
the total mass. Results for 1 million particle simulations with low (left) and high (right) resolution grids are displayed. The dashed line 
is the Kepler scaling of orbits, r 3 oc (tg — t) 2 . At constant internal mass, each orbit converges asymptotically to Kepler scaling under 
adequate resolution. 



5.1 Oblate and prolate spheroids 

To construct spheroidal distributions, we squeezed or 
stretched the axes of spheres to achieve the sought geometry. 
We consider two cases, a prolate r = spheroid of initial 
axes = (2, 1, 1) and an oblate one with axes = (2, 2, 1). We 
then performed iV-body calculations with 10 4 , 10 5 and 10 6 
particles and compared the outcome with the pressure-less 
fluid solutions for the same initial configurations. 

The results are displayed on Fig. |(| This graphs the 
gravitational energy W as a function of time. The dimen- 
sionless axes xt of the fluid solution are also displayed, where 
we have indicated the formation of discs or spindles, accord- 
ing to whether a single or two axes vanished at the time 
indicated. Spindle or disc each gives rise to sharp features in 
the profile of W . The numerical solutions for r = or 1 show 
clear dependencies on particle numbers, in the sense that the 
larger- N calculations map the features of the fluid solution 
more closely. Better agreement with the fluid solutions are 
to be expected as we increase particle number, since the fluid 
solutions corresponds to N — > oo. This is difficult to assess 
quantitatively for the solutions as a whole. However we may 
isolate features that support this view. For instance, as N is 
increased from 10 4 to 10 6 particles, the maximum potential 
energy achieved in both cases displayed increases and cor- 
responds to 'spindles' in the fluid solution. For the case of 
the t = 1 (oblate) fluid spheroid, a spindle forms following 
collapse of the major-axis at t ~ 0.623. The time of maxi- 
mum potential is t = 0.646, 0.642 and 0.638 respectively for 
the N = 10 4 , 10 s and 10 6 runs, when these maxima shifts 
upward with N, from 9.57 to 13.2, and 15.7 for the largest- 
N run. The value for the fluid solution — > oo formally. Note 
that the N = 10 6 run is the only one with a rapid re-collapse 
to a disc singularity (see inset at t m 0.64) similar to the fluid 
solution. A similar comment can be made for the r = (pro- 
late) fluid spheroid, where the sequence of singularities is 



reversed: a spindle forms first, followed by disc and spindle 
singularities. The two numerical calculations both follow the 
fluid solution relatively well, with one important difference: 
at the time of the first spindle, t ~ 0.370, the 10 4 particle 
run shows an increase in W/W(0) much reduced compared 
with the 10 6 particle run (3.77 to 5.27, or 70% as much); 
at t ~ 0.525 a second spindle forms but now the two N- 
body runs develop very similar extrema in W/W(0): 5.42 
and 5.61, respectively. This means that following rebound 
through the firsts spindle, both calculations suffer a compa- 
rable degree of orbit mixing, which then smoothes out the 
second singularity seen in the fluid. As for the r = 1 case, 
the phase of the iV-body curves tunes up to the fluid solu- 
tion as N increases: thus the first singularity at t ~ 0.370 
is found at t = 0.385 (N = 10 4 ) and t = 0.376 (N = 10 6 ), 
which differs with the fluid solution by only 1.6%. 

The two cases r = and 1 displayed on Fig. ^| both 
develop spindles which will reach arbitrarily large potential 
energy as N — > oo. Both behave in a qualitatively identical 
way in this respect. We decided therefore to investigate in 
more details only the relation of the solution for the oblate 
t = 1 case to the number of particles, before exploring tri- 
axial initial configurations. 



5.2 Series of oblate spheroids 

We constructed oblate spheroidal distributions as indicated 
above. The equator of the spheroids lies in the XY plane. 
The diameter was kept fixed and only one axis resized to 
achieve the desired aspect ratio, hence the gravitational en- 
ergy W is magnified with decreasing a3(0)/oi(0). In the 
limit 0,3(0) -^Owe compute a radial free-fall time 

t«(a 3 [0] =0) = J^x t a (Eq. |lj) ~ 0.5116 . 
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Figure 5. (A) Collapse factor C versus particle number N. Uniform spheres follow the scaling relation oc N 1 / 3 (solid line) well when 
the grid mesh resolves particles individually (black circles). The open circles are results obtained by ALP+88 with a direct-summation 
code. Crosses were taken from Theis & Spurzem (1999) for Plummer models using the direct-integration code NBODY6++. Their data 
follows the same power-law as uniform distributions (dotted line). 



Figure 5. (B) Collapse factor C versus particle number N for uniform spheroidal distributions. The collapse factor for axisymmctric 
spheroids is well fitted with the empirical power JV 1 / 6 (solid line). The A^/^-scaling law of spheres is also shown for reference, for 
both uniform- and Plummer- initial distributions (dotted lines). The open triangles are results of Boily et al. (1999) obtained with a 
direct-summation ./V-body algorithm. The horizontal arrow indicates a collapse factor »s O(20), such that 50% of all triaxial ellipsoids of 
infinite particle number (i.e. pressure-less fluid) collapse to smaller values. 



The time of maximum contraction would therefore lie be- 
tween this and the value (|2l]) for spheres. 

The results are listed in Table |^. The errors on the col- 
lapse factor C — W/W(0) are computed from variations 
about the mean value for fixed number of particles. The 
aspect ratios initially lie between 1/6 and 9/10, for particle 
numbers ranging from 10 4 to 10 7 . All results are graphed 
as triangles on Fig. |B|(b). (The effect of different initial as- 
pect ratios on C is discussed in §5.3.) We have added points 
obtained from simulations with direct-summation codes by 
Boily et al. (1999) for small-A systems to those of the 
present study (filled triangles on the figure). For large- A'' 
calculations (N > 10 5 and beyond) the black triangles mark 



a gently increasing trend, well matched with a power-law 
dependence C oc N a with a « 1/6. For N w 10 4 or smaller, 
the fit remains good but note the large scatter for points 
obtained with a direct-summation code (open triangles on 
the figure). The range of collapse factors measured for the 
10,000 particle runs listed in Table ^ allows for a multiplica- 
tive factor of 3/2 between maximum and minimum values of 
W/W(0). The results for the direct-summation runs would 
allow a somewhat larger range, of perhaps 5/2. Whether this 
is cause for concern is debatable because of the small num- 
ber of runs in this bin; it may be that particle-particle in- 
teractions, better resolved in the direct-summation scheme, 
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Figure 6. Evolution of the self-gravitating energy versus time for prolate (r = 0, left) and oblate (r = 1, right) spheroids. The solid lines 
give the solution for a pressure-less homogeneous fluid; broken lines are TV-body realisations. The labels give the grid cell and particle 
numbers. On the top panels, the set of curves shown at the bottom are the dimensionless axes xi integrated from (^) andjl6|) and we 
have indicated the times when spindles or discs form. Note how cusps develop in the analytic solution each time a disc or spindle forms. 



Table 2. Collapse factor C = W m ax /W (0) for uniform spherical distributions. The models have varying particle number (N), mesh size 
and linear resolution (m and () in model units, but the same initial total potential energy Vl^O) = 1.20. The free-fall time % = t ma x 
corresponds to the time when W reached a maximum; the analytic value (fell) is given in round brackets. At that time Lagrangian radii 
enclosing 30% and 60% of the mass were measured; their values are given here respective to their initial values, L o (30%) = 0.670 and 
L o (60%) = 0.843. 

N m I iff = tmax L /L Comments 

W (0) 

~1Q- 3 (0.7854) 30% 60% 



10 4 


64 


1.6 


0.795 


29.5 


35.1 


31.9 




10 4 


64 


1.6 


0.799 


29.6 


35.1 


31.9 


Reduced St 


10 4 


32 


3.1 


0.798 


18.6 


26.1 


18.9 


Lower resolution 


10 s 


128 


0.8 


0.787 


67.0 


85.9 


75.0 




10 s 


64 


1.6 


0.788 


67.6 


91.8 


89.2 




10 5 


64 


1.6 


0.788 


56.1 


73.8 


64.4 


St X 2 


10 s 


32 


2.6 


0.786 


32.9 


58.9 


41.0 


Lower resolution 


10 6 


128 


0.8 


0.784 


97.3 


133.8 


112.4 




10 6 


64 


1.6 


0.781 


73.4 


121.6 


89.7 


Lower resolution 


; x io 7 


128 


0.8 


0.784 


225.9 


352.31 


290.0 





increase the scatter somewhat, though not the mean values, 
which we recover well with the FFT scheme. 

The trend with particle number TV for spheroidal dis- 
tributions is never well fitted with the scaling oc N 1 ^ 3 of 
spherical distributions, though the data differ by only small 
factors for small- N systems. The results of Section [] suggest 
that any increasing trend of collapse factor C with particle 
number would be a sensitive function of the symmetry of the 
initial distribution, or of evolution towards axisymmetry in 
the course of evolution. 



5.3 Another look at the LMS flow 

We investigated the role played by the Lin-Mestel-Shu insta- 
bility in numerical iV-body calculations of violent relaxation. 
For homogenous systems, the LMS instability develops as an 
aspherical system collapses from rest first down its shortest 
axis. For oblate spheroids, the collapse down the minor axis 
z occurs in a time 



*(* = «)« J SteCEq© , (22) 



12 Boily, Athanassoula & Kroupa 



Table 3. Collapse factor C = W m ax/W(0) for uniform spheroidal distributions. Symbols as for Table | 

N m I t m ax — 7-7 7777^ Comments 

ai(0) W(0) 
-MO" 3 < 1 



10 4 


32 


3.1 


0.737 


4:5 


12.27 


10 4 


32 


3.1 


0.732 


4:5 


11.04 


If) 4 


32 


3.1 


0.697 


2:3 


10.23 


If) 4 


64 


3.1 


0.695 


2:3 


12.40 


If) 4 


32 


3.1 


0.659 


1:2 


11.52 


10 4 


64 


3.1 


0.650 


1:2 


12.65 


10 4 


32 


3.1 


0.646 


1:2 


9.23 


10 4 


32 


3.1 


0.608 


1:3 


11.88 


10 4 


64 


1.6 


0.602 


1:3 


17.26 


10 4 


32 


3.1 


0.574 


1:6 


8.85 


10 4 


32 


3.1 


0.584 


1:6 


10.50 


10 s 


32 


3.1 


0.722 


4:5 


17.44 


10° 


64 


1.6 


0.724 


4:5 


17.86 


10 5 


128 


0.8 


0.563 


4:5 


18.09 


10 s 


32 


3.1 


0.686 


2:3 


13.32 


10 r ' 


32 


3.1 


0.642 


1:2 


13.29 


10 s 


32 


3.1 


0.597 


1:3 


18.40 


10 5 


64 


3.1 


0.595 


1:3 


24.1 


10 5 


128 


1.6 


0.597 


1:3 


30.1 


10° 


32 


3.1 


0.558 


1:6 


14.54 


10 s 


64 


3.1 


0.559 


1:6 


15.69 


10 s 


64 


1.6 


0.563 


1:6 


17.19 


10 5 


128 


0.8 


0.563 


1:6 


17.57 


10 8 


128 


0.8 


0.752 


9:10 


33.17 


10 6 


128 


0.8 


0.682 


2:3 


17.68 


10 6 


64 


0.8 


0.634 


1:2 


15.60 


10 8 


128 


0.8 


0.637 


1:2 


15.84 


10 (i 


128 


0.8 


0.595 


1:3 


44.16 


10 8 


128 


0.8 


0.593 


1:3 


47.42 


10 8 


128 


O.i 


0.593 


1:3 


47.48 


10" 


64 


1.6 


0.596 


1:3 


28.89 


10 8 


128 


0.8 


0.556 


1:6 


23.20 


10 7 


128 


0.8 


0.718 


4:5 


29.33 


10 7 


128 


0.8 


0.680 


2:3 


22.55 


10 7 


128 


0.8 


0.595 


1:3 


57.45 


10 7 


128 


0.8 


0.553 


1:6 


36.92 



<W/W >= 11.6 ±1.3 



<W/W >= 18 ± 3 



Higher resolution 
Lower resolution 
<W/W >= 28 ± 11 



<W/W >= 37 ± 10 



where tg is the free-fall time for spheres. For the initially 
flattest spheroids in our sample, of as(0) : ai(0) = 1 : 6, we 
compute t = 0.41 tg ~ 0.32. The data given in Table |^ and 
Fig. g]& ^show that maximum collapse (or, W) occurs at 
later times, as collapse down the major axis sets in. This 
applies to all simulations. Thus self-induced LMS-type of 
instabilities are not sufficient by themselves to halt collapse 
since the kinetic energy dispersion a grows anisotropically. 
This anisotropy persists up to the time of maximum W and 
beyond, even for low-iV particle number calculations, and 
imprints the virialised equilibrium that follows (Boily et al. 
1999). 

A question remains which concerns the relative impor- 
tance of the LMS instability compared to the fragmentation 
modes of instability that control collapse in spherical 
symmetry (cf. Section 2.1). Both types of instability will 
develop during the collapse of a spheroid, however the 
growth rate of the LMS instability deduced from (J22I) will 



be higher for spheroids with initially small aspect ratios. In 
Section 2.2, we have argued that the growth rate of the LMS 
instability is always more rapid than the fragmentation 
instability if the initial aspect ratio is less than about 0.75. 
We would, therefore, expect a signature of this instability 
in the form of a stream at a later stage of collapse, namely 
the bounce. We seek out evidence for this in our sample of 
runs of Table |§|. Below we refer to the stream as an 'LMS 
flow', which should not be confused with the instability 
described by Lin and co-workers. 

Since ( |22[ ) shows a relation between minor-axis collapse 
time and initial aspect ratio, we expect a similar relation 
between maximum collapse factor C = W max /W(0) (the 
bounce) and initial aspect ratio for the cases when an LMS 
flow drives the dynamics at that time. Crucial to the ar- 
gument is the relative phase of the velocity components at 
the bounce. For the chosen spheroidal initial conditions, the 
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solution from which the histogram was constructed does not 
suffer from any type of instability, the poor agreement with 
the data would suggest that instabilities other than the LMS 
instability are not completely negligible to set the system 
properties at the bounce. 

We note that for TV = 10 4 particle runs and initial as- 
pect ratio > 1/2 or so, the data lie near the normalised C = 1 
value (cf. Table | and Fig. |, open squares), hence are not 
caught in the sinusoidal pattern driven by the streaming 
motion. However our data also indicates that 10 4 particle 
runs with initial aspect ratios < 1/2 or so match the range 
of values obtained with larger- TV runs. Thus fragmentation 
modes of instability may still be the dominant factor con- 
trolling the bounce for simulations with particle numbers 
< 10 4 and initially near-spherical distributions, but not so 
when the initial aspect ratio is sufficiently small. 



Figure 7. Normalised collapse factor C versus initial aspect ratio 
for all runs. The values C taken from Table |i| were normalised to 
the mean of their respective TV series. The open squares are the 
results for N = 10 4 ; black triangles for all others. The dashed line 

is the function sin (9 + 9 — 2) , where 9 = 2iz — . The histogram 

ai 

at the bottom shows the relative phase of the minor-axis and 
major-axis velocities around the time when a singularity forms. 

motion can be divided in cylindrical z and R components, 
and we may align the minor axis with the z-component of 
the reference frame. The expectation for LMS flows is that 
when both <v z > and <vr> are negative inward (in phase) 
at the bounce, the value of W achieved should be larger than 
when <v z > and <vr> are of opposite signs (out of phase), 
i.e. one inward, the other outward. Unless fragmentation or 
other types of instability manage to erase the signature LMS 
flow, the relative phase of the velocity components at the 
bounce will be set by the initial system aspect ratio. 

We may compare the ensemble of calculations of Ta- 
ble gj, first by normalising individual values of W/W(0) for 
given TV to the mean value for that series; in this way we 
remove the scaling oc TV 1 / 6 between simulations with differ- 
ent particle numbers. We then sort the normalised values by 
increasing order of the initial aspect ratio, as(0)/or(0). The 
results are displayed on Fig. The sinusoidal pattern of the 
data is unmistakable. The data may be fitted with an sinu- 
soid of amplitude m 1/2, which is comparable or larger than 
the intrinsic scatter of the points at given aspect ratio. Thus 
the LMS flow has equal or more impact on the potential and 
the system configuration than other factors which predict 
no dependence with initial aspect ratio, such as the growth 
of velocity dispersion by internal fragmentation modes (see 
Section 2). At the bottom of Fig. |^ we have sketched the 
relative phase of the velocity components obtained for the 
analytic fluid solution for oblate spheroids. The step-wise 
histogram indicates in-phase (high step) motion or out-of- 
phase (low step) motion. The arrows indicate the polarity 
of the motion. The correspondence with the numerical data 
is only suggestive: the pattern itself appears somewhat out 
of phase with the histograms. Since the pressure-less fluid 



6 RESULTS FOR TV-BODY TRIAXIAL 
COLLAPSE (t/OOR 1) 

We now extend our study to triaxial configurations. We re- 
peated the exercise of Section |j| but this time varying two 
initial aspect ratios as parameters. The set-ups used to ob- 
tain numerical and analytic solutions were as before. The 
results are illustrated on Fig. ^. We considered three triaxial 
configurations with r = 1/3, 1/2 and 2/3, so covering both 
prolate and oblate structures. For these cases the analytic 
fluid solution did not develop spindles and hence the poten- 
tial energy remained finite for the duration of integration. 
As for the spheroidal calculations, the numerical TV-body 
calculations come ever closer to the analytic solution with 
increasing particle number. For example, for TV = 10 6 parti- 
cle runs, the potential energy already comes within 20% of 
the fluid solution at maxima. Higher particle numbers would 
only bring modest differences and convergence as TV — > oo is 
therefore very slow. The 10 4 particles runs remains approxi- 
mately 50% out of step with the fluid values, and hence the 
quantities involved with such low- TV calculations of relax- 
ation processes are to be treated with caution in applications 
to galaxy or halo formation problems. 



7 DISCUSSION AND CONCLUSION 

We have sought to constrain the tidal field developing 
around galaxies and haloes as they form. To do this we stud- 
ied the growth of gravitational energy during the violent 
relaxation of ellipsoidal bodies. We used both an analytic 
pressure-less gas model and TV-body numerical integration 
to set absolute limits on calculations of galaxy and halo for- 
mation involving TV point masses. 

We found using the pressure-less gas model that close to 
9/10 (86%) of all ellipsoidal triaxial configurations increase 
their gravitational energy by at most a factor 40 (cf. Figs, [j] 
and ^). We confirmed this with TV-body calculations using 
up to one million particles. 

We studied axisymmetric and spherical uniform distri- 
butions. We extended the scaling of collapse factor with 
particle number for homogeneous spheres, C oc TV 1//3 , to 
TV = 16 millions with the code SUPERBOX. We noted that 
axisymmetric spheroidal distributions also show increasing 
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Figure 8. Evolution of the self-gravitating energy W 
versus time for three ellipsoidal configurations of initial 
morphology given by r (cq. Jis|]). In all cases, the solid 
line represents the pressure-free fluid solution, while 
broken lines are the TV-body realisations with parti- 
cle numbers as indicated. The evolution of the dimen- 
sionless ellipsoid shown at the bottom for 

reference. The TV-body calculations come ever closer to 
the fluid solution with increasing particle number. Note 
that the fluid solution reaches a finite maximum W in 
all cases displayed. 



collapse factor with particle number: The run of data points 
is well fitted by the power-law C oc TV 1 / 6 , much gentler than 
for spherically symmetric distributions (see Fig. ^[b]). The 
scatter in the data for spheroids as such does not allow to 
fall back on the C oc N 1 ^ 3 scaling for spheres, even when we 
try matching runs with different initial aspect ratio (cf. Ta- 
ble ^) . We pointed out that the extrema of binding energy 
and hence tidal forces met by such systems depend sensi- 
tively on the formation of prolate structures (spindles) and 
are much gentler otherwise. 

The growth of velocity dispersion during collapse can 
be attributed both to global fragmentation modes and 
Lin-Mestel-Shu-type of pancaking. We presented evidence 
to the effect that in calculations involving more than 
N — 10 4 particles, the Poissonian seeds of fragmentation 
modes leads to growth in kinetic energy such that the sum 
remains smaller than, though not negligible before, the 
growth of kinetic energy attributable to the surface mode 
(Lin-Mestel-Shu) of instability. For N < 10 4 and initial 



aspect ratio > 1/2, the data suggests that fragmentation 
modes play an equally important role for up to N — 10 4 
particles (see Fig. |). For N > 10 s , the LMS flow at 
the bounce sets the system properties both in terms 
of potential depth and velocity components for the full 
range of the initial conditions studied here, where aspect 
ratios where taken in the range from 1/6 to 9/10. For 
this range of particle number and more, the collapse of 
systems with different initial morphologies can be distin- 
guished without ambiguity on Fig. 5, which gives confidence 
on that the gravatitational infall has been properly resolved. 

Taken together, these results imply that the forma- 
tion of axially- or spherically-symmetric haloes and galax- 
ies lead to deeper potentials during the violent relaxation 
phase and hence to more pronounced tidal fields than for 
non-symmetric ones. The stronger tidal fields would in turn 
reduce the rate of survival of sub-condensations or satel- 
lites orbiting within them. Consequently, we would expect 
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galactic halo morphology and satellite populations to be 
correlated, in the sense that galaxies with rounder massive 
haloes would harbour fewer satellites, while triaxial haloes 
would harbour many more satellites, all other parameters 
being equal. This does not take into account the long-term 
fate of galactic satellites: Tidal forces do not subside after 
virialisation, and eventually will cause the disruption of all 
galactic satellites after a period of time (e.g. Ibata et al. 
1994, Klessen & Kroupa 1998, Bullock et al. 2001). Thus 
the above statement refers to the time of formation only, 
once the system has virialised. 

Another direct consequence of our results is that spher- 
ically symmetric haloes should be more centrally concen- 
trated than non-spherical ones in virial equilibrium. We may 
expect this to bear on the kinematics of observed galaxies. 
Unfortunately the current observational constraints on the 
halo shapes are not sufficiently precise to allow us to test our 
prediction. Halo axial ratios have been measured for hardly 
over a dozen galaxies. What is more worrisome, however, is 
that the different techniques seem to give systematically dif- 
ferent results. Merrifield (2002) summarises nicely the situa- 
tion for disc galaxies. Their halos appear to be axisymmetric 
and oblate, with their axes of symmetry coalligned with the 
disk axes. The most reliable measurements for the minor to 
major axis ratios come from polar rings, but galaxies hav- 
ing such structures may not be a representative sample of 
disk galaxies, because they might be the results of recent 
mergers. Measurements from the flaring of the HI disk give 
systematically smaller values than those of polar rings mea- 
surements. An application to our own Galaxy (Oiling and 
Merrifield 2000) shows that such values may be valid only if 
the value of the distance from the sun to the Galactic center 
and the local Galactic rotation speed are smaller than what 
is currently believed. Certainly some progress is necessary 
before the measurements attain the precision we need for 
testing our prediction. 

Gravitational lensing offers some hope by constraining 
the distribution of total (dark + visible) gravitational mass 
inside a given volume from the symmetry of the lensed 
image, which will not respect the centre of mass of the 
system if it is not spherically symmetric. For instance, 
Mailer et al. (2000) have applied such a lensing technique 
to the spiral B1600+434. The deconvolution procedure 
however does suggest that the shape of the halo deduced 
remains dependent on the choice of halo density profile 
(isothermal or otherwise) and symmetry. It may be that 
the systematic application of such techniques to sufficiently 
large samples would reveal a correlation in the sense that 
we indicated above. 

The results obtained for uniform-density distributions 
should be contrasted with results obtained for non-uniform 
initial distributions. Theis & Spurzem (1999) investigated 
the morphological evolution of initially cold Plummer 
distributions. The collapse factors they obtained are given 
on Fig. |E|(a) and |E](b). These confirm earlier results by 
ALP+88 of lower values of C for non-uniform systems 
and our own arguments of Section 2. The curve fitting 
the Theis & Spurzem data is shown on Fig. ^|(b) shifted 
down with respect to the one obtained for uniform spheres. 
This new curve now intersects with the collapse factors 
obtained for aspherical distributions and N ~ 10 3 particles. 



A Plummer model shows an extended envelope of mass 
density p oc r _5//2 . Thus for the same particle number, 
the collapse factor of a Plummer model is reduced in 
comparison with a uniform sphere, presumably due to shell 
crossing taking place near the centre. However, Plummer 
models with larger particle numbers also collapse by 
larger factors, and hence other systems with initially steep 
profiles will, too. In spherical symmetry, the collapse of 
mass distribution with radial dependence (e.g. p oc r~ a ) 
has bearing on accretion problems, since the mass shells 
reach the centre a various rates in time. The currently 
favoured road to galaxy formation would have many clumps 
converging to the centre of mass. Since the scaling we have 
obtained for uniform-density profiles may be extended to 
non-uniform profiles, as shown with the Plummer model, 
we may hope that the relation of gravitational gradient to 
initial morphology will also find application to cosmological 
models of galaxy and galactic halo formation. 

To recover the physics of collapsing systems adequately 
in TV-body calculations requires a sufficiently large number 
of particles in order to disentangle effects of mass distribu- 
tion and morphology. The results of Fig. ^(b) suggest a fidu- 
cial number TV ~ 100, 000 particles as a clean demarcation 
(where a large gap appears between the lower dotted and 
solid line on the figure) between initially spherical, axisym- 
metric and triaxial distributions. Furthermore, the effect of 
varying the initial density profile appears only to shift the 
zero-point of the curves, and hence does not affect the rela- 
tion of maximum collapse factor C = W ma x/W(0) to particle 
number. 
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